C.................... RSSUBC.FOR ................................
C::::::::::::::::::::::::::: Module FIELD_LINE_APEX_PARAMS :::::::::::::::::
C.... Contains parameters associated with the Apex field line grid
C.... Written by P. Richards June-July 2021.
      MODULE FIELD_LINE_APEX_PARAMS !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
        IMPLICIT NONE
        !.. These are the related to the global solution grid. 
        !.. glatX, glonX, and altX are reset by Apex routines
        !** altmax has little influence on the grid setup cpu
        REAL altmaX,glatX,glonX,altX,hrX,SIGNLAT
        Data  altmaX, glatX,  glonX ,  altX  ,   hrX     !.. PGR
     >      /99000.0, 43.0, -71.5 , 250.0 , 00.0   /
      END MODULE FIELD_LINE_APEX_PARAMS
C::::::::::::::::::::::::::: Module FIELD_LINE_APEX_GRID :::::::::::::::::
C.... Contains Apex parameters associated with the FLIP field line grid
C.... Written by P. Richards June-July 2021.
      !.. qdlat,qdlon       ! quasi-dipole latitude and longitude (deg)
      !.. bmagF             !.. magnetic field strength
      !.. gdlat,gdlon       ! geodetic latitude and longitude (deg)
      !.. DIPX,DECX,COSDIPX !  Magnetic dip, declination and cos(dip)
      MODULE FIELD_LINE_APEX_GRID
        IMPLICIT NONE
        INTEGER IDIM
        parameter (IDIM=401)  !.. Check dimensions are the same in all modules
        real gdlat(IDIM),gdlon(IDIM),bmagF(IDIM),qdlat(IDIM),qdlon(IDIM)
        REAL DIPX(IDIM),DECX(IDIM),COSDIPX(IDIM)
      END MODULE FIELD_LINE_APEX_GRID
C::::::::::::::::::::::::::::::::::::::: RFIELD ::::::::::::::::::::::::::::
C....... this program is responsible for calling the subroutine that sets up
C....... the dipole magnetic field grid and filling the field parameter arrays,
C....... altitude, gravity etc. SeeG.J. Bailey, Planet. Space Sci., Vol. 31,
C....... No. 4, pp. 389-409, 1983.
      SUBROUTINE RFIELD(IWR,PCO,RE,Z0,SCAL,SCALK,ZLBDY,Q,VPERP,
     &             VTOT,COSDIP,JB,JR,JNQ,JTI,ZHI,ZLO,JVERT,JVERC)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER JB,JR,JNQ,JTI,JVERT,JVERC,JVOL
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,RE,GLATD(401),GLOND(401)
      DIMENSION FD(9),COSDIP(401), X(401),Q(401)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      common/ExB/eq_vperp,div_vem(401)

      RPTS=(JMAX+1)/2-1         ! just converting to the real
      DH=1.0/RPTS               ! distance between adjacent points in x coord.
      ALTR=200.0+200.0*(F107+F107A)/2.0/70. ! used in limiting flux calculation
      eq_vperp=vperp

      JVOL=1
      DO 10 J=1,JMAX
        CALL FIELD(J,JMAX,PCO,RE,SCAL,SCALK,Q(J), X(J),FD)
        JU=JMAX+1-J                ! symetric point of J about equator plane
        IF(J.EQ.1) SREF=FD(8)*1.E5 ! distance along magnetic field
        Z(J)=FD(1)                 ! altitude in km
        IF(Z(J).LE.ZLBDY)  JB =J   ! lower boundary
        IF(Z(J).LT.ALTR)   JR =J   ! index for limiting flux
        IF(Z(J).LT.ZHI) JNQ=J      ! index for plasmaspheric fluxes
        IF(Z(J).LT.120.) JVOL=J      ! index for plasmaspheric fluxes
        GL(J)=FD(4)                ! latitude
        BM(J)=FD(6)                ! magnetic field strength
        BG(J)=FD(7)                ! coordinate factor
        COSDIP(J)=FD(3)            ! cos(dip)
      !..     IF(JTI.GT.1) GO TO 10
        GR(2,J)=-FD(5)             ! gravity
        R(2,J)=(RE+FD(1))*1.E5     ! radial distance to grid point
        IF(J.NE.1) SL(J)=SREF-FD(8)*1.E5   ! distance to grid point along B
      !..  Evaluate the ExB term for the continuity equation
        CALL GET_DIVDRIFT_VELOCITY(RE,PCO,X(J),VPERP, DIV_VEM(J))
        IF(DABS(GL(J)).LT.1.0D-4)     GO TO 20
      !..  Set up field line parameters in Southern hemisphere
        GL(JU)=-GL(J)
        Z(JU)=Z(J)
        BM(JU)=BM(J)
        BG(JU)=BG(J)
        COSDIP(JU)=-COSDIP(J)
        GR(2,JU)=-GR(2,J)
        SL(JU)=SREF+FD(8)*1.E5
        R(2,JU)=R(2,J)
        DIV_VEM(JU) = DIV_VEM(J)
  10  CONTINUE
  20  CONTINUE

      JQ=(JMAX+1)/2
      IF(Z(JQ).LT.99) THEN
        WRITE(6,696) NINT(Z(JQ)) 
 696    FORMAT(/'   ** ERROR: Field line apex altitude is',I3,' km.'
     >      ,1X,'It does not reach the ionosphere ***'/)
        CALL RUN_ERROR    !.. print ERROR warning in output files
        STOP 
      ENDIF

      IF(ZLO.LT.Z0.OR.Z0.GE.120.) ZLBDY=Z0

       !--- Find altitude points for setting lat, long for vertical grid
      JVERT=1
      DO 245 IL=1,JMAX/2
        IF(Z(IL).LT.250) JVERT=IL
 245  CONTINUE

      JVERC=JMAX+1-JVERT   !... conjugate index

      !.. Determine Apex grid values on the FLIP field line grid
      CALL FLIP_APEX_GRID(JMAX,Z)
      !.. Transfer Apex grid values to the FLIP grid
      CALL UPDATE_FLIP_GRID(JMAX,Z,COSDIP,BM,GL,GLATD,GLOND,GR)
      
      !... Output the flux tube parameters if requested
      IF(IWR.NE.7) RETURN   !.. No output if not needed in FLIP PRINT
      REWIND (17)  !... overwrite previous values if flux tube moves

      !.. Calculate total flux tube volume for header
      VTOT=0.0
      DO J=JVOL+1,JMAX
        VTOT=VTOT+BM(JVOL)*(SL(J)-SL(J-1))/BM(J)
      ENDDO

      !-- Write the flux tube parameter headers
      WRITE(17,100) IDAY/1000,MOD(IDAY,1000),F107,F107A,AP(1)
     >   ,SEC/3600.0,PCO,BLON,NINT(SL(JMAX)*1.0E-5),VTOT

      TOTVOL=0.0
      DELS=0.0
      DELZ=0.0
      VELEM=0.0
      DO J=1,JMAX
        IF(J.GT.1) THEN
          DELS=(SL(J)-SL(J-1))/1.0E5
          DELZ=ABS(Z(J)-Z(J-1))
          IF(J.GT.JVOL) THEN
            VELEM=BM(JVOL)*(SL(J)-SL(J-1))/BM(J)
            TOTVOL=TOTVOL+VELEM
          ENDIF
        ENDIF
        WRITE(17,110) J, Z(J), (SL(J)*1.0E-5), DELZ ,DELS,
     >   (ACOS(COSDIP(J))*57.3),COSDIP(J),BM(J),TOTVOL,GR(2,J),
     >   (R(2,J)*1.0E-5),DIV_VEM(J),GL(J)*57.296,GLATD(J),GLOND(J),VELEM
      ENDDO
 110  FORMAT(I5,5F10.2,2F9.2,1P,E10.2,0P,2F9.1,1P,E10.2,0P,3F7.2,
     >  1P,E10.2)

      RETURN
 100  FORMAT(/9X,' * * * * File for field line grid parameters,
     > EUV fluxes cross sections * * * * *'
     > ,//,'  YEAR= ',I4,'    DAY= ',I3,'   F107= ',F5.1,'   F107A= '
     >   ,F5.1,'     Ap= ',F5.1
     > ,//9X,' **** The field line grid parameters at UT=',F8.2,' ****'
     > ,//,'  L-shell =',F7.3,';     Magnetic longitude= ',F5.1
     > ,/,' Length of flux tube = ',I8,' km'
     > ,'     Total flux tube volume = ',1P,E10.2,' cm-3'
     > ,//,'    pt     alt    arc_len   alt_dif   arc_dif   dip_ang'
     > ,2X,'cos_dip  mag_fld  Tube_vol  gravity   Re+alt  Div_Vem'
     > ,3X,'B_lat  G_lat G_long  DEL_VOL')
      END
C:::::::::::::::::::: GET_DIVDRIFT_VELOCITY ::::::::::::::::::::::::::::::
      SUBROUTINE GET_DIVDRIFT_VELOCITY (RE,PCO,X,VPERP, DIV_VEM)
C
C.... This subroutine calculates the ExB term in the continuty equation.
C.... This expression comes from Bailey, PSS 1983, page 392. Note that Re in
C.... his equation (5) is the equatorial radius to the flux tube
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      REAL RE
 
      NUMERATOR = 6.0000 * ((DSIN(X))**2.0)*(1.0+(DCOS(X))**2.0)
      DENOMINATOR = RE*PCO*1.0E5*((1.0+3.0*(DCOS(X))**2.0))**2.0
      DIV_VEM = NUMERATOR/DENOMINATOR
      RETURN
      END
 
C:::::::::::::::::::::::::: QFIELD ::::::::::::::::::::::::::::::::::::::::
        SUBROUTINE QFIELD(JMAX,PCO,RE,Z0,SCAL, SCALK,Q)
C
C.... This program determines Q values.
C.... JMAX = # of grid points, PCO = lshell, Z0 = lower boundary,
C.... scal = scaling factor to ensure that the x-coord ranged from 0 to 1.
C.... SCALK & Q are returned.
C.... Reference: G.J. Bailey, Planet. Space Sci., Vol. 31, No. 4,
C               pp. 389-409, 1983.
C
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      DIMENSION Q(401)
      REAL RE
 
      !... R0 = equatorial radius to flux tube.
      R0=RE*PCO
      PTS=(JMAX+1)/2
      IPTS=PTS
      RAT=(RE+Z0)/RE
      QMAX=(SQRT(1-RAT/PCO))/(RAT**2)
      SCALK=1/SINH(SCAL*(QMAX))
      DO J = 1, JMAX
        DX=1-(J-1)/(PTS-1)
      !..  q=the dipole coord  determined from -- sinh(kq)=dx
        SCX=DX/SCALK
        Q(J) = DLOG(SCX+DSQRT(SCX**2+1))/SCAL
      ENDDO
 
      RETURN
      END
 
C:::::::::::::::::::::::::: FIELD ::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE FIELD(J,JMAX,PCO,RE,SCAL,SCALK,Q, XR,FD)
C..........................................................
C  this program determines the grid point spacing given just jmax=
C   # of grid points, pco=lshell, z0=lower boundary, scal=scaling
C   factor;; x=colatitude, the field parameters are transferred
C   through fd(i);; r0=equatorial radius to flux tube, initial values
C   for x and r are set. dh=distance between points in the x coordinate
C  Reference: G.J. Bailey, Planet. Space Sci., Vol. 31, No. 4,
C               pp. 389-409, 1983.
C--------------------------------------------------
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
C     temp_store is used to store previous computed values
      DIMENSION FD(9)
      REAL RE
      X=.5
      R0=RE*PCO
      !... Values for equatorial point
      IF(J .EQ. (JMAX+1)/2) THEN
         X=1.570796327
         SHX=1
         CHX=0
      ELSE
         JITER=0

         !... loop for Newton solver
  10     CONTINUE
         SHX=DSIN(X)
         CHX=DCOS(X)
      !.. the next 6 lines are a newton solver for the equation f(x)=0
      !.. this determines the colatitude x
         FEX=(R0**2)*(SHX**4)-(RE**2)*CHX/Q
         DEX=(R0**2)*4*(SHX**3)*CHX+(RE**2)*SHX/Q
         X=X-FEX/DEX
         JITER=JITER+1
         IF(JITER.GT.99) THEN
           WRITE(6,*) 
     >      ' ** ERROR: Terminate after 100 iterations by Newton solver'
           CALL RUN_ERROR    !.. print ERROR warning in output files
           STOP
         ENDIF
         IF(ABS(FEX/(DEX*X)).GT.1.0E-6) GO TO 10
      ENDIF
      !..end newton solver

      !... XR is a returned value other than X in order to keep previous value of
      !... X, generated during last CALL, from being modified by current CALL with
      !... some random value.
      XR = X
      !..  r=radial distance to grid point. sqth=numerator of mag field
      R=R0*SHX**2
      SQTH=DSQRT(3*(CHX**2)+1)
      !..  fd(1)=altitude. fd(2)=sin(dip). fd(3)=cos(dip). fd(4)=latitude
      FD(1)=R-RE
      FD(2)=2*CHX/SQTH
      FD(3)=SHX/SQTH
      FD(4)=1.570796327-X
      !..  fd(5)=gravity. fd(6)=mag field strength. fd(7)= coord factor.
      FD(5)=FD(2)*3.98E+10/((RE+FD(1))**2)
      FD(6)=8.271E+25*SQTH/(R*1.E+5)**3
      FD(7)=SCAL*DCOSH(SCAL*Q)*SCALK*1.E-5*RE**2*SQTH/R**3
      !..  ds & fd(8)= arc length along field calc by 2 methods
      XX=DLOG(1.732*CHX+SQTH)
      FD(8)=R0*.28868*(XX+DSINH(XX)*DCOSH(XX))
 
      RETURN
      END 
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE BTOG(IEXP,BLOND,BLATD,GLATD,GLOND)
C....... Conversion from centered dipole to geographic coordinates. (PLAT
C....... ,PLON)= coords of north magnetic pole. The geomagnetic longitude
C....... is measured east from the geomagnetic prime meridian at 291 degrees
C....... east geographic.
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !....... degrees to radians
      BLONR=BLOND/57.29578
      BLATR=BLATD/57.29578
      !...... Position of point in geomagnetic coords
      XM=COS(BLATR)*COS(BLONR)
      YM=COS(BLATR)*SIN(BLONR)
      ZM=SIN(BLATR)
      !...... rotate coords in north-south direction
      XG=XM*SIN(PLAT)+ZM*COS(PLAT)
      YG=YM
      ZG=-XM*COS(PLAT)+ZM*SIN(PLAT)
C
      !..... geographic latitude and longitude converted to degrees
      GLATD=ASIN(ZG)*57.29578
      GLOND=(PLON+ATAN2(YG,XG))*57.29578
      IF(GLOND.GE.360) GLOND=GLOND-360
      RETURN
      END
C::::::::::::::::::::::::: EPHEM ::::::::::::::::::::::::::::::::::::::
      SUBROUTINE EPHEM(IDAY,ETRAN)
      DIMENSION EPTRAN(26)
      DATA EPTRAN/4338,4374,4398,4404,4392,
     >  4374,4344,4320,4302,4296,4302
     > ,4320,4338,4356,4356,4344,4320
     > ,4290,4260,4236,4218,4224,4248,4284
     > ,4320,4350/
      ANUM=(MOD(IDAY,1000)+15.)/15.
      INUM=ANUM
      FRAC=ANUM-INUM
      ETRAN=(EPTRAN(INUM)*(1.0-FRAC)+EPTRAN(INUM+1)*FRAC)*10.
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SYM(JMAX,HE,EHT,UN,Z,SZA)
C........... This routine gives symmetric ambient parameters
C....... NOTE WELL ! This program does not produce symmetry at night
C....... because the density at grazing incidence is not symmetric
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      DIMENSION UN(401),Z(401),HE(401),EHT(3,401),SZA(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      DO 25 J=1,JMAX/2
      JU=JMAX+1-J
      ON(JU)=ON(J)
      O2N(JU)=O2N(J)
      N2N(JU)=N2N(J)
      HN(JU)=HN(J)
      HE(JU)=HE(J)
      PHION(JU)=PHION(J)
      SZA(JU)=SZA(J)
      TN(JU)=TN(J)
      UN(JU)=-UN(J)
      DO 25 I=1,4
      IF(I.LE.3) TI(I,JU)=TI(I,J)
      IF(I.LE.3) EHT(I,JU)=EHT(I,J)
      IF(I.GT.0) N(I,JU)=N(I,J)
 25   CONTINUE
      UN(JMAX/2+1)=0.0D0
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE BRACE(H,NI,TE,TB,RAT)
C...... Empirical Te model from Brace and Theis GRL p275, 1978
      DOUBLE PRECISION H,NI,TE,TB,RAT
      IF(H.LT.120.OR.H.GT.500) RETURN
      TB=1051+(17.07*H-2746)
     > *EXP(-5.122E-4*H+6.094E-6*NI-3.353E-8*H*NI)
      RAT=TE/TB
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SETDT(DT,DTSET,DTMAX,JTIDT,ITAU,Z,JMAX,ITF,ITER,TWIDT)
C......... routine for setting time step DT ....
      DOUBLE PRECISION DT,DTSET,DTMAX,Z(401),TWIDT,GL,DTSAV
      DIMENSION D(19),T(2)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      DATA JTICNT /0.0/, TWISTRT/1.8/
      !...... negative TWIDT changes time to begin/end twilight time step
      IF(TWIDT.LT.0) THEN
         TWISTRT= 2.0
         TWIDT = -TWIDT
      ENDIF
      
      !------- Output diagnostics to log file
      IF(ITER.LT.-2) WRITE(6,626) ITER
 626     FORMAT(2X,'*** Difficulty converging, ITER=',I3
     >  ,' - time step reduced')

      DTADD=INT(0.5*DT)
      !...      IF(DT.GT.300.AND.DTADD.LT.300) DTADD=300.
      IF(ITER.LE.-2) GO TO 161
      JTICNT=JTICNT+1
      IF(ITER.LT.6.AND.JTICNT.GT.JTIDT) DT=DT+DTADD
      IF(DT.GT.DTMAX) DT=DTMAX
      GO TO 163
C
      !.............. adjust DT for NON-convergence .......
 161  ITAU=ITAU-DT
      SEC=SEC-SNGL(DT)
      DT=NINT(0.5*DT)
      JTICNT=0
      IF(JTIDT.LT.3) JTIDT=3

      !....... advance program time step and ut
 163   IHDT=(INT(DT)/30)*30
       IF(DT.GT.60.) DT=FLOAT(IHDT)
       IHDT=(INT(DT)/60)*60
       IF(DT.GT.120.) DT=FLOAT(IHDT)
       IHDT=(INT(DT)/300)*300
       IF(DT.GT.600.) DT=FLOAT(IHDT)
       IF(DT.GT.DTMAX) DT=DTMAX
       ITAU=ITAU+DT
       IF(ITF.EQ.1.OR.ITF.LT.0) SEC=SEC+DT
C
      !........ determine local time at both ends of flux tube .......
       CALL AMBS(1,Z(1),GL(1),D,T,CHIN,SATN,GLONN,GLATN)
       CALL AMBS(JMAX,Z(JMAX),GL(JMAX),D,T,CHIF,SATF,GLONF,GLATF)
      !....... ensure reasonable time step at sunrise and sunset ....
      IF(ITF.NE.1.AND.ITF.GE.0) RETURN
      IF(DT.LE.TWIDT) RETURN
      !...      IF(ABS(SATN-12).LT.4.AND.ABS(SATF-12).LT.4) RETURN
      IF(CHIN.LE.1.5.AND.CHIF.LE.1.5) RETURN
      IF(CHIN.GE.TWISTRT.AND.CHIF.GE.TWISTRT) RETURN
      IF(CHIN.GE.TWISTRT.AND.CHIF.LE.1.5) RETURN
      IF(CHIF.GE.TWISTRT.AND.CHIN.LE.1.5) RETURN
C
      !........... reduce time step to 1800 seconds near sunrise and sunset
      DTSAV=DT
      IF(DT.GT.1800) DT=1800
      IF(DT.GT.DTMAX) DT=DTMAX
      ITAU=ITAU-DTSAV+DT
      SEC=SEC-SNGL(DTSAV)+SNGL(DT)
C
      !............ set time step to TWIDT at sunset
       CALL AMBS(1,Z(1),GL(1),D,T,CHIN,SATN,GLONN,GLATN)
       CALL AMBS(JMAX,Z(JMAX),GL(JMAX),D,T,CHIF,SATF,GLONF,GLATF)
      IF(CHIN.LE.1.5.AND.CHIF.LE.1.5) RETURN
      IF(CHIN.GE.TWISTRT.AND.CHIF.GE.TWISTRT) RETURN
      DTSAV=DT
      IF(DT.GT.TWIDT) DT=TWIDT
      IF(DT.GT.DTMAX) DT=DTMAX
      ITAU=ITAU-DTSAV+DT
      SEC=SEC-SNGL(DTSAV)+SNGL(DT)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PHEAD(IFILE,PCO,GLONN,GLONF,GLATN,GLATF,ALT)
C....... Subroutine for printing basic information in each of the
C....... output files
      DOUBLE PRECISION PCO,GL,ALT
      INTEGER FEXISTS
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/SOL/UVFAC(59),EUV
      CHARACTER*60 HEMIS
      DATA FEXISTS/-9/
      
            WRITE(IFILE,676)
 676  FORMAT(/2X,'^^^ FLIP MODEL:- 2023-10-13 version ****'
     > ,1X,'Check the bottom of *LOG*.txt file for the latest updates.'
     > ,//2X,'*** Consult either Phil Richards (pgrichds@gmail.com)'
     > ,1X,'or Pat Dandenault (Patrick.Dandenault@jhuapl.edu)'
     > ,1X,'for the appropriate level of recognition for model use.')

      IF(IFILE.EQ.3) THEN
         HEMIS = 'the NORTHERN hemisphere * * * * * * * * * *'
      ELSE IF(IFILE.EQ.4) THEN
         HEMIS = 'the EQUATORIAL PLANE * * * * * * * * * *'
      ELSE IF(IFILE.EQ.9) THEN
         HEMIS = 'the SOUTHERN hemisphere * * * * * * * * * *'
      ELSE IF(IFILE.EQ.6) THEN
         HEMIS = 'BOTH hemispheres * * * * * * * * * *'
      ELSE
         HEMIS = ' * * * * * * * * * *'
      ENDIF
      WRITE(IFILE,112) HEMIS
 112  FORMAT(//5X,'**** Initial Parameters in the ',A60)
 114  FORMAT(/1X,' YEAR=',I4,',  DAY=',I3,',  L-SHELL=',F8.3,
     >  ',  F10.7=' ,F6.1,',  F10.7A=',F6.1,',  AP=',F6.1)
      WRITE(IFILE,114) IDAY/1000,MOD(IDAY,1000),PCO,F107,F107A,AP(1)
      
      !... print lats and longs as integers or else reals
      IF(IDAY.GT.999) THEN
         WRITE(IFILE,115) NINT(GLATN),NINT(GLONN),-NINT(GLATF)
     >  ,NINT(GLONF),NINT(ALT)
      ELSE
         WRITE(IFILE,915) GLATN,GLONN,-GLATF,GLONF,ALT
      ENDIF
 115  FORMAT(/'  Geographic (lat,long):- North (',I3,'N,',I5,'E)'
     > ,2X,',  South (',I3,'S,',I5,'E)  @',I6,'km altitude')
 915  FORMAT(/'  Geographic (lat,long):- North (',F5.1,' N,',F6.1,' E)'
     > ,',  South (',F5.1,' S,',F6.1,' E)  @',F9.2,' km altitude')

      IF(FEXISTS.EQ.-9) CALL TFILE(42,FEXISTS) !.. see if fluxes from file
      IF(FEXISTS.NE.0) THEN
         WRITE(IFILE,87) 
 87      FORMAT(/3X,'+++ EUV fluxes read from file +++'/)
      ELSEIF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3) THEN
        WRITE(IFILE,116) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 116    FORMAT(/'   Using HEUVAC EUV. Initial scaling factors: 303.78='
     >  ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >  ', 1025.72=',F5.1)
      ELSEIF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4) THEN
        WRITE(IFILE,117) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 117    FORMAT(/' Using TIMED-SEE EUV. Initial scaling factors: 303.78='
     >   ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >   ', 1025.72=',F5.1)
      ELSEIF(NINT(UVFAC(58)).GE.0) THEN
        WRITE(IFILE,118) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 118    FORMAT(/'   User provided EUV scaling factors: 303.78='
     >  ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >  ', 1025.72=',F5.1)
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PROFIN(ISOL,J,JMAX,TALT,F107,N,TI,TN,SL,GR,ON,HN,HPEQ
     > ,HEPRAT,XIONN,IVLPS)
C....... set up rough initial o+, h+, and temperature profiles
      DOUBLE PRECISION ALT,N(4,401),TI(3,401),TN(401),SL(401),GR(2,401)
     > ,ON(401),HN(401),HPEQ,XHOLD,HEPRAT,XIONN(9,401),TALT
      IF(ISOL.NE.0) GO TO 10
      ALT=TALT
      IF(ALT.LT.50) ALT=50
      !..... Temperatures approximated by a log function
      TI(3,J)=904.0*DLOG(ALT)-3329.0
      IF(TI(3,J).GT.5000) TI(3,J)=5000
      TI(1,J)=TI(3,J)
      TI(2,J)= TI(3,J)
      IF(J.GT.(JMAX+1)/2) GO TO 10
      !...... Insert initial O+ profile: first low altitudes
      XHOLD=-5.5E-4*(ALT-250.0)**2
      IF(XHOLD.GT.70.0)XHOLD=70.0
      IF(ALT.LE.500) N(1,J)=(F107/74.0)*5E5*EXP(XHOLD)
      !....... now high altitudes
      XHOLD=(SL(J+1)-SL(J))*1.9E-7*GR(2,J)/(TI(3,J)+TI(1,J))
      IF(XHOLD.GT.70.0)XHOLD=70.0
      IF(ALT.GE.275) N(1,J)=N(1,J-1)*EXP(XHOLD)
      N(1,JMAX+1-J)=N(1,J)
C------- N+ density is 10% of O+
      IF(IABS(IVLPS).EQ.11) THEN
         XIONN(1,J) = 0.1 * N(1,J)
         XIONN(1,JMAX+1-J)= 0.1 * N(1,J)
      ELSE
         XIONN(1,J) = 0.0
         XIONN(1,JMAX+1-J)= 0.0
      ENDIF
      IF(HPEQ.EQ.1001) WRITE(6,90) NINT(ALT),TI(3,J),N(1,J)
 90   FORMAT(1X,'PROFIN',I7,F8.0,1P,E10.2)
      !...... Initial H+ and He+ profiles
 10   IF(HPEQ.LT.0) RETURN
      IF(ON(J).GT.0) N(2,J)=N(1,J)*HN(J)/ON(J)
      IF(N(2,J).GT.HPEQ.OR.ALT.GE.2000) N(2,J)=HPEQ
      N(4,J)=0.0
      IF(IABS(IVLPS).GE.9) N(4,J)=HEPRAT*N(2,J)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE DATRD1(ISOL,ITF,ITFS,I1,I2,IPP,MAXIT,JTIMAX,DTSET
     > ,DTMAX,JTIDT,TWIDT,ZLBDY,IREAD,HPEQ,FPAS,HFAC,IVLPS,TSTOP,JWR
     > ,ZLO,ZHI,IWIND,IVIBN2,IODDN,INPLUS,IHEP,HEPRAT,ITEMP,IHPOP
     > ,IHPKP,BLATPR)
      INTEGER ISCH
      DOUBLE PRECISION DTSET,DTMAX,TWIDT,ZLBDY,HPEQ,FPAS,HFAC,TSTOP
     > ,ZLO,ZHI,HEPRAT,DHDU_FAC,COLFAC,CFAC,RTS(199),BLATPR
      COMMON/CFACT/DHDU_FAC,COLFAC,CFAC(9),ISCH(9)
      SAVE RTS
      
      !.. TSTOP no longer checked to be > 38 hours. December 2008
      READ(5,*) TSTOP     !.. Length of FLIP run in hours. 
      TSTOP=ABS(TSTOP)    !.. kept for safety

C------ Test to see if file exists.
      CALL TFILE(15,IOPEN)
      IF(IOPEN.EQ.0) THEN
         WRITE(6,79) 
 79      FORMAT(/1X,' ** ERROR: FLIPRUN.DDD file does not exist ***'/)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF
C------ Read from file until a good value for ITEMP - get rid of headers
 7    READ(15,*,ERR=7) ITEMP
         IF(ITEMP.LT.0.OR.ITEMP.GT.1) ITEMP=1
         I1=ITEMP
      READ(15,*) IHPOP
         IF(IHPOP.LT.0.OR.IHPOP.GT.1) IHPOP=1
         I2=I1+IHPOP
         IF(ITEMP.EQ.0.AND.IHPOP.EQ.1) I1=2
         IF(ITEMP.EQ.0.AND.IHPOP.EQ.1) I2=2
      READ(15,*) IHEP
         IF(IHEP.LT.0) IHEP=1
      READ(15,*) INPLUS
         IF(INPLUS.LT.0) INPLUS=1
      READ(15,*) IODDN
         IF(IODDN.LT.0) IODDN=1
      READ(15,*) IVIBN2
         IF(IVIBN2.LT.0) IVIBN2=1
      !.. Switch to enhanced O+ + N2 reaction rate if no N2v
      RTS(199)=0
      IF(IVIBN2.EQ.0) RTS(199)=-1
      CALL RATS(1,999,999,999,RTS)

      READ(15,*) TDHDU         !.. To alter the Wind DHDU factor
      DHDU_FAC=0.1        
      IF(TDHDU.GT.-.0001.AND.TDHDU.LT.0.0001) THEN
        DHDU_FAC=0.1        
      ELSEIF(TDHDU.GE.0.05.AND.TDHDU.LE.2.0) THEN
        DHDU_FAC=TDHDU
        WRITE(6,'(A,F8.3)') '  *** CAUTION DH/DU was changed to',TDHDU
      ELSEIF(TDHDU.LT.0.05.OR.TDHDU.GT.2.0) THEN
         IF(TDHDU.LT.0.05)
     >   WRITE(6,'(A,F9.3)') ' ** ERROR: DHDU is < 0.05, = ',TDHDU
         IF(TDHDU.GT.2.0)
     >   WRITE(6,'(A,F8.3)') ' ** ERROR: DHDU is > 2.00, = ',TDHDU
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF

      READ(15,*) BLATPR  !.. Latitude for special print 2020-08-30

      !....... Now set IVLPS to determine which species to solve
         IVLPS=0
         IF(IODDN.GE.1) IVLPS=1
         IF(IVIBN2.GE.1) IVLPS=2
         IF(IHEP.NE.0.AND.IVLPS.EQ.0) IVLPS=-9
         IF(IHEP.NE.0.AND.IVLPS.GT.0) IVLPS=+9
         IF(INPLUS.NE.0.AND.IVLPS.LE.0) IVLPS=-11
         IF(INPLUS.NE.0.AND.IVLPS.GT.0) IVLPS=+11

      READ(15,*) DTSET
         IF(DTSET.LT.1.OR.DTSET.GT.1800) DTSET=1800
      READ(15,*) DTMAX
         IF(DTMAX.LT.1.OR.DTMAX.GT.7200) DTMAX=1800
         IF(DTSET.GT.DTMAX) DTSET=DTMAX
      READ(15,*) TWIDT
         IF(TWIDT.LT.10.OR.TWIDT.GT.3600) TWIDT=600
         IF(TWIDT.LT.300) WRITE(6,109) 
 109     FORMAT(/' **** WARNING: TWIDT < 300, Pyyddd file may be'
     >   ,' very large'/)
      READ(15,*) IWIND
      READ(15,*) HPEQ
        IHPKP=0
        IF(NINT(HPEQ).EQ.-1) IHPKP=1
        IF(HPEQ.GT.1E4) HPEQ=-HPEQ
c        IF(ISOL.NE.0.AND.HPEQ.GE.0.0) HPEQ=-HPEQ
      READ(15,*) HEPRAT
        IF(HEPRAT.LT.0.001) HEPRAT=0.001
        IF(HEPRAT.GT.2.0) HEPRAT=2.0
      READ(15,*) FPAS
        IF(FPAS.LT.0.0) WRITE(6,90)
 90     FORMAT(//2X,'FPAS <0, plasmaspheric heating prop. to TEC')
        IF(FPAS.GT.1.0) WRITE(6,91)
 91     FORMAT(//2X,'FPAS >0, no conjugate e* & no plasmaspheric heat')
      READ(15,*) HFAC
        IF(HFAC.LT.0) WRITE(6,92) ABS(HFAC)
 92     FORMAT(//'  HFAC<0, multiplying ionospheric heating rate by +'
     >   ,F5.1)
      !...
      READ(15,*) JWR
         IF(JWR.LT.0) JWR=0
      READ(15,*) IPP
         IF(IPP.LT.0) IPP=0
         IF(JWR.NE.0.AND.IPP.EQ.0) IPP=1
      READ(15,*) ZLO
         IF(ZLO.LT.70.0) ZLO=300.0
      READ(15,*) ZHI
         IF(ZHI.LT.70.0) ZHI=300.0
      READ(15,*) ITF
      !....         IF(ITF.LT.1) ITF=1
         ITFS=ITF
         IF(ISOL.EQ.0.AND.ITF.GT.0) ITF=3
      READ(15,*) JTIMAX
         IF(JTIMAX.LT.0) JTIMAX=9999999
         !...IF(ISOL.EQ.0.AND.JTIMAX.LT.333) JTIMAX=333
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C..... This subroutine reads in the geophysical parameters for the run
      SUBROUTINE DATRD2(ISOL,IWIND,AP,BLON,F107,F107A,IDAY,SEC,UVFAC)
      REAL AP(7),UVFAC(59),ZFLUX(37),DLAT,DLON
      DATA RUNHRS/0.0/, DLAT/99.0/

      READ(5,*) AP(1)
         IF(ABS(AP(1)).GT.400) THEN
            WRITE(6,*) ' ** ERROR: AP is too big = ',AP(1)
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF

      READ(5,*) F107
         IF(F107.LT.60.OR.F107.GT.500) THEN
             WRITE(6,*) ' ** ERROR: F107 is out of range =',F107
             CALL RUN_ERROR    !.. print ERROR warning in output files
             STOP
         ENDIF
      READ(5,*) F107A
         IF(F107A.LT.60.OR.F107A.GT.500) THEN
             WRITE(6,*) ' ** ERROR: F107A is out of range =',F107A
             CALL RUN_ERROR    !.. print ERROR warning in output files
             STOP
         ENDIF
      READ(5,*) JDAY
         IF(ISOL.EQ.0.OR.JDAY.GT.0) IDAY=IABS(JDAY)
         IF(IDAY.GT.0.AND.IDAY.LT.367) IDAY=IDAY+2000000
         IF(IDAY.LT.1900000) IDAY=IDAY+1900000
         IF(IDAY.LT.1000000.OR.IDAY.GT.9999999) THEN
           WRITE(6,*) ' ** ERROR: Day number out of range = ',IDAY
           CALL RUN_ERROR    !.. print ERROR warning in output files
           STOP
         ENDIF

      !.. Determine if the year is in the FLIP magnetic pole range. 
      !.. Negative DLAT triggers the test
      CALL MPOLE_INC(IDAY,-DLAT,DLON) 

      !...... The universal time (SEC) may be changed in DATRD3
      READ(5,*) STRTHRS
         IF(ISOL.EQ.0)  SEC=STRTHRS*3600

      DO I=1,59
        UVFAC(I)=1.0
      ENDDO
      !....... Read in the UV factors
      READ(15,*) (UVFAC(I),I=1,37)
      READ(15,*) (UVFAC(I),I=38,50)

      !.. Multiplication factors for neutral densities
      !.. Species order = 51-He, 52-O, 53-N2, 54-O2, 55-Ar, 57-H, 58-N
      !.. 56 = total mass density. 59 = Hot O density
      READ(15,*) (UVFAC(I),I=51,59)

      !.. See if UV factors need to be adjusted. UVFAC(58) is used to
      !.. switch between different fluxes
      UVFAC(58)=0.0
      IF(NINT(UVFAC(1)).GT.0.0001) THEN
         UVFAC(58)=2    !.. user input EUV factors
      ELSEIF(NINT(UVFAC(1)).EQ.-2) THEN
         UVFAC(58)=-2   !.. Use SEE factors
         CALL FACSEE(F107,F107A,UVFAC,ZFLUX)
      ELSE
         UVFAC(58)=-1   !.. Use HEUVAC factors
         CALL FACEUV(F107,F107A,UVFAC,ZFLUX)
      ENDIF

      !.. See if user input S-R factors
      IF(UVFAC(38).LE.0.0) THEN 
        UVFAC(58)=UVFAC(58)-2    !.. default S-R factors
        CALL FACSR(UVFAC,F107,F107A)
      ENDIF

      !... check the magnitude of the euv fluxes ....
      IUV=0
      DO 67 I=1,45
 67   IF(UVFAC(I).GT.99.OR.UVFAC(I).LT.0) IUV=1
      IF(IUV.EQ.1) THEN
        WRITE(6,*)
     >  '  **** EUV or S-Runge factors may be wrong - check FLIPRUN.DDD'
        WRITE(6,751) (UVFAC(I),I=1,50)
 751    FORMAT(10F6.1)
      ENDIF

      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE DATRD3(ISOL,BLON,SEC,DTSET,F107,JMAX,PCO,Z0
     > ,SCAL,GRD,L_STAG,L_SHORT,L_PHASE,IDAY)
      use apex     !.. Richmond's Apex Magnetic Field Model
      INTEGER IDAY
      DOUBLE PRECISION PCO,Z0,SCAL,DTSET,XHOLD
      REAL L_STAG,L_SHORT,L_PHASE,GLOND,GLATD,GLONP
      REAL LVALUE,BLONG
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !...... If not a new field line then read dummy variables so they
      !...... are not altered
      IF(ISOL.NE.0) THEN
         READ(15,*) DUMVAR
         READ(5,*)  L_SHORT
         READ(5,*)  DUMVAR
         READ(5,*)  L_STAG
         READ(5,*)  L_PHASE
         READ(15,*)  DUMVAR
         READ(15,*)  DUMVAR
         READ(15,*)  DUMVAR
         RETURN
      ENDIF
 
      !....... else setting up a new field line, input parameters
      !...... JMAX is the number of grid points on field line
      READ(15,*) JMAX
      !..... PCO (P coordinate) is the L-SHELL or geographic latitude. 
      READ(5,*) PCO
      !..... BLON is lonitude - magnetic OR geographic if PCO is latitude
      READ(5,*) BLON

      !.. stagnation point and phase for tear drop model
      READ(5,*) L_STAG
      IF(NINT(L_STAG).NE.0) THEN
         WRITE(6,'(//4X,A,/4X,A,A,/4X,A/)') 
     >  ' *@$%* @$%* @$%* @$%* @$%* *@$%* @$%* @$%* @$%* @$%* @$%*',
     >  ' ***** This Apex version of FLIP does not yet support',
     >  ' ExB drifts *****',
     >  ' +++++ Contact pgrichds@gmail.com if you need ExB drifts ++++'
         STOP
      ENDIF
      READ(5,*) L_PHASE

      !.. If PCO large it must be geographic latitude need to convert
      !.. to L-shell and magnetic longitude. Modified 2022-08-15
      IF(PCO.LT.1.02.OR.PCO.GT.9) THEN
         GLATD=PCO
         GLOND=BLON

         !.. Determine if the highest altitude on field line is in the ionosphere.
         !.. If so, it doesn't make sense to use 250 km as the field line altitude
         !.. The 250 km puts the field line near hmF2 above a station
         CALL APEX_LVALUE(IDAY,GLATD,GLOND,250.0,LVALUE,BLONG)
         
         !.. Set up the Apex global grid and get L-value. 
         IF((LVALUE-1.0)*RE.GE.600.0) THEN
            CALL APEX_LVALUE(IDAY,GLATD,GLOND,250.0,LVALUE,BLONG)
         ELSE
            CALL APEX_LVALUE(IDAY,GLATD,GLOND,120.0,LVALUE,BLONG)
         ENDIF
         PCO=LVALUE     !.. Update L-value to Apex/IGRF value

      ELSE  !.. Else treat PCO and BLON as L-shell and magnetic longitude 
         GLONP=BLON-74.5
         IF(GLONP.LT.0) GLONP=GLONP+360
         !.. locate magnetic pole
         CALL MPOLE_LOC(GLONP,50.0,IDAY,PLAT,PLON)
         !...    transform magnetic to geographic coordinates---
         BLAT=57.2958*ACOS(SQRT(1.0/REAL(PCO)))
         !.. Use PCO as L-value to determine magnetic latitude and 
         !.. set up Apex Coordinates
         CALL LVTOGX(IDAY,100.0,REAL(PCO),BLON)
      ENDIF

      IF((BLON-55.0)**2.0+(GLATD+30.0)**2.0.LT.25.0*25.0)
     >  WRITE(6,'(4X,A,/4X,A,/4X,A)') 
     >  'WARNING: This flux Tube is in the South Atlantic Anomaly.',
     >  'The displaced dipole may not be a good representation.',
     >  'Not a serious problem for the local ionosphere.'
     
      !.. L_SHORT is where the trajectory cuts the short axis
      IF(L_STAG.GT.1.15) L_SHORT=PCO

      !.. Give a warning for N2v at low latitudes
      IF(NINT(L_STAG).GT.0) THEN
        IF(PCO.LT.1.15.OR.L_SHORT.LT.1.15)
     >  WRITE(6,'(/4X,A,/4X,A)') 
     >  'WARNING: For this flux Tube, L gets less than 1.15.',
     >  ' Turn off N2 vibration in FLIPRUN.DDD (IVIBN2=0)'
      ENDIF

      !-- correct JMAX if illegal value and make sure it is an odd number
      IF(JMAX.LT.150) JMAX=55*PCO+125
      JMAX= 2 * (JMAX/2) + 1
      IF(JMAX.LT.299) JMAX=299
      IF(JMAX.GT.399) JMAX=399
      IF(NINT(L_STAG).NE.0) JMAX=399   ! max pts for ExB drift

      !.. Z0 is the absolute lower boundary altitude
      READ(15,*) Z0
         !.. reset Z0 for ExB drift
         IF(INT(L_STAG).EQ.0) THEN
            IF(Z0.LT.80) WRITE(6,*) ' *** Z0 reset to 80 km'
            IF(Z0.LT.80) Z0=80
         ELSE
            Z0=90
            WRITE(6,97) NINT(Z0)
 97         FORMAT('  *** Z0 reset to',I4,' km for ExB drift')
         ENDIF
         IF(Z0.GT.190.0) THEN
            Z0=190.0
            WRITE(6,*) ' Z0 too high - reset to 190 km'
         ENDIF
C
      !....... SCAL distributes the grid points along B
      READ(15,*) SCAL
         IF(SCAL.LT.0.5.OR.SCAL.GT.11) SCAL=(8.0/PCO+2.0)
         IF(NINT(L_STAG).NE.0) SCAL=3.5   !.. SCAL for ExB drift
      !........ switch for variable (-1) or constant altitude steps in vertical
      READ(15,*) GRD
      !..   IF(GRD.LT.0) GRD=-1
C
      !...... Set up default starting time at noon by finding the approx. UT
      !...... IF  ABS(SEC) < 24 assume it is LT in hours to start  
       IF(SEC.LT.0) CALL SET_NOON(SEC,PCO,BLON)
C
 95    FORMAT(I6,9F10.3)
      !...... Default initial time step for new field line
      XHOLD=40*(250/F107)**1.5
      IF(DTSET.GT.XHOLD) DTSET=INT(40*(250/F107)**1.5)

      RETURN
      END

C::::::::::::::::::::::::::::: SET_HPEQ ::::::::::::::::::::::::::::::::
C... This function sets the default H+ density at the equator. This is done
C... if the current HPEQ is less than 30 (can be negative)
      SUBROUTINE SET_HPEQ(PCO,HPEQ)
      IMPLICIT NONE
      DOUBLE PRECISION PCO,HPEQ,HP_FULL
      !... Take user value or set default H+ density at equator.
      IF(HPEQ.LT.30) THEN
         !... set default HPEQ
         HP_FULL=1.22E5/PCO**4
         !... if value between 0.1 and 4, take fraction of HP_FULL
         IF(HPEQ.LT.0.1.OR.HPEQ.GT.4.0) HPEQ=1.0
         HPEQ=DABS(HPEQ)*HP_FULL
         IF(HPEQ.LT.1) HPEQ=1
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::: DATRD4 ::::::::::::::::::::::::::::::::
C........ Used for reading additional run control parameters
      SUBROUTINE DATRD4(ISOL,TIMEON,TIMEOF,AURFAC)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL TIMEON,TIMEOF,AURFAC,AMP_RING,RC_TIME_CON,TMAX,SIGT
     >   ,LMAX,SIGL,PROP
      !....... COLFAC= O+ - O collision frequency factor. CFAC, and ISCH
      !....... are for future expansion
      COMMON/CFACT/DHDU_FAC,COLFAC,CFAC(9),ISCH(9)
      COMMON/RC/AMP_RING,TMAX,SIGT,LMAX,SIGL,RC_TIME_CON,PROP
      !... common DD transfers EDDY to CALODDN and CALVN2
      COMMON/DD/D(7,VD),E(7,VD),B1(7,VD),B2(7,VD),RDZ(VD),EDDY

      !...... This is the O+ - O collision factor
      READ(15,*) COLFAC
         IF(COLFAC.LT.0) COLFAC=1.0
         IF(COLFAC.LT.0.5.OR.COLFAC.GT.4.0) THEN
            WRITE(6,*) ' ** ERROR: COLFAC out of order = ',COLFAC
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF

      READ(15,*) EDDY
         IF(EDDY.GT.1.0E7) THEN
            WRITE(6,*) ' ** ERROR: EDDY diff. coefficient too big',EDDY
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF

      !.. Amplitude of the ring current heating
      READ(15,*) AMP_RING
      AMP_RING=AMP_RING*5.0E-05
      IF(AMP_RING.LE.0.0) AMP_RING=0.0

      !.. time of maximum and sigmas for ring current heating
      READ(15,*) TMAX,SIGT
      READ(15,*) LMAX,SIGL

      !... proportion of ring current heating going to electrons
      READ(15,*) PROP
      IF(PROP.GT.1.0) PROP=1.0
      IF(PROP.LT.0.0) PROP=0.0

      READ(15,*) RC_TIME_CON     !.. decay time for ring current (hrs)

      !.. input on time and off time for aurora and energy factor
      READ(15,*) TIMEON,TIMEOF
      !.. 2020-11-06: allow extra heating to straddle midnight
      IF(TIMEON.GE.0.0.AND.TIMEON.LT.24.0.AND.TIMEOF.LE.TIMEON.AND.
     >    TIMEOF.LT.24.0) TIMEOF=TIMEOF+24.0
      READ(15,*) AURFAC
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE CENFOR(IEXP,BLOND,BLATD,ALT,CF)
C....... Calculation of centrifugal force. Note that it is asymmetrical with
C....... respect to the geomagnetic equator
      DOUBLE PRECISION ALT,BLATD,CF,RP,CL,DIP,CHI
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !....... Find geographic latitude and longitude
      CALL BTOG(IEXP,BLOND,SNGL(BLATD),GLATD,GLOND)
      RP=(RE+ALT)*1.0E5
      !..... rotate the centrifugal force vector in longitude
      XG=COS(GLATD/57.296)*COS(GLOND/57.296-PLON)
      YG=COS(GLATD/57.296)*SIN(GLOND/57.296-PLON)
      ZG=0.0
      !...... rotate in latitude
      XM=XG*SIN(PLAT)-ZG*COS(PLAT)
      YM=YG
      ZM=XG*COS(PLAT)+ZG*SIN(PLAT)
C
      CL=(90.0-BLATD)/57.296
      DIP=ASIN(2*DCOS(CL)/DSQRT(3.0*DCOS(CL)**2+1.0))
      CHI=CL-DIP
      !....... Centrifugal force =CF
      CF=RP*ANGVEL*(DCOS(CHI)*(COS(BLOND/57.296)*XM+SIN(BLOND/
     > 57.296)*YM)-DSIN(CHI)*ZM)
      RETURN
      END
C:::::::::::::::::::::::::::::::: FACEUV :::::::::::::::::::::::
      SUBROUTINE FACEUV(F107,F107A,UVFAC,ZFLUX)
C----- This routine uses the EUV scaling from Richards et al.[1994]
C----- The EUVAC flux model is based on the F74113 solar reference
C----- spectrum and Hinteregger's scaling factors. This subroutine
C----- just provides the scaling factors as a function of the proxy
C----- (F107+F107A)/2. Modified in March 2009 to return fluxes as 
C-----  well as factors.
      IMPLICIT NONE
      INTEGER I,L,LMAX
      REAL F107,F107A,F107AV,A,B
      REAL UVFAC(59),HFG200(37),ZFLUX(37),ZFX(37)
      !..Ratio of EUVAC fluxes for P=200 and P=80 (day 1974113)
      DATA HFG200/2.202,1.855,2.605,3.334,1.333,17.522,4.176,4.0
     >  ,1.4,3.694,1.791,5.385,1.889,1.899,3.427,2.051,1.392,1.619
     >  ,1.439,2.941,1.399,2.416,1.512,1.365,1.570,1.462,2.537,1.393
     >  ,1.572,1.578,1.681,1.598,1.473,1.530,1.622,1.634,1.525/
      DATA LMAX/37/
      !--  Test to see if need to scale - see DATRD2 subroutine      
      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3) THEN
         !........... EUV scaling
         F107AV=(F107+F107A)*0.5
         DO I=1,LMAX
           A=(HFG200(I)-1)/120.0
           B=1-A*80.0
           UVFAC(I)=A*F107AV+B
           IF(UVFAC(I).LT.0.8) UVFAC(I)=0.8
         ENDDO
      ENDIF

      !.. These fluxes, which are divided by 1.0E+9, are for day 1974113.
      !.. Note!!!! fluxes doubled below 250A, while the shortest wavelength 
      !.. fluxes have been tripled
      DATA ZFX/2.4665,2.1,3.5,1.4746,4.4,3.,3.537,1.625,.758,.702,
     > .26,.17,.141,.36,.23,.342,1.59,.53,.357,1.27,.72,.452,.285,
     > .29,.383,.314,.65,.965,6.9,.8,1.679,.21,.46,3.1,4.8,.45,1.2/

      !.. Determine the actual EUV fluxes. Note that for historical reasons, 
      !.. the UV factors and fluxes run in opposite directions
      DO L=1,LMAX
         ZFLUX(L)=ZFX(L)*1.0E+9*UVFAC(LMAX+1-L)
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::::: FACSEE :::::::::::::::::::::::
C----- This routine produces scaling factors for the SEE solar EUV 
C----- fluxes. It uses P=(F107+F107A)/2. The scaling is based on a solar  
C----- minimum spectrum for day 2002092 (P=196) and day 2007284 (P=68)
C----- Written by P. Richards January 2009
      SUBROUTINE FACSEE(F107,  !.. Daily F10.7 cm flux
     >                 F107A,  !.. 81-day average F10.7 cm flux
     >                 UVFAC,  !.. IN/OUT Scaling factors for EUV flux
     >                 ZFLUX)  !.. OUTPUT: Solar EUV fluxes
      IMPLICIT NONE
      INTEGER I,L,LMAX
      REAL F107,F107A,F107AV,A,B
      REAL UVFAC(59),SEEFAC(37),ZFLUX(37),ZFX(37)
      !.. ratio of fluxes for day 2002092 (P=196) and day 2007284 (P=68)
      DATA SEEFAC/1.82,1.50,1.51,1.78,1.57,14.39,1.84,3.71,1.43,8.02,
     >  1.64,2.94,1.96,1.73,1.90,2.19,1.40,2.15,1.03,3.02,1.51,2.15,
     >  1.65,1.30,1.78,1.27,2.24,1.35,1.29,1.72,1.93,1.85,1.38,2.48,
     >  1.50,1.53,2.58/
      DATA LMAX/37/
      !--  Test to see if need to scale - see DATRD2 subroutine      
      IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4) THEN
         !........... EUV scaling
         F107AV=(F107+F107A)*0.5
         DO I=1,LMAX
           A=(SEEFAC(I)-1)/128.0
           B=1-A*68.0
           UVFAC(I)=A*F107AV+B
           IF(UVFAC(I).LT.0.8) UVFAC(I)=0.8
         ENDDO
      ENDIF

      !.. These are SEE fluxes for day 2007284 (P=68), divided by 1.0E+9
      DATA ZFX/3.106,2.1,3.5,1.888,4.4,4.433,4.734,2.42,0.731,0.7,0.26,
     >  0.1832,0.353,0.36,0.5454,0.712,0.795,0.265,0.349,0.791,0.5,
     >  0.775,0.817,0.29,0.649,1.051,0.65,0.4824,6.3736,0.739,4.51,
     >  0.21,0.46,3.558,6.26,0.28,0.4721/

      !.. Determine the actual EUV fluxes. Note that for historical reasons, 
      !.. the UV factors and fluxes run in opposite directions
      DO L=1,LMAX
         ZFLUX(L)=ZFX(L)*1.0E+9*UVFAC(LMAX+1-L)
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::: FACSR ::::::::::::::::::::::::::::::::
C.... The Schumann-Runge factors are scaled according to F10.7
C.... Using SEE fluxes
C.... P. Richards March 2011
      SUBROUTINE FACSR(UVFAC,F107,F107A)
      DIMENSION UVFAC(59),SRFLUX(8),SEEFAC(8)
      !.. ratio of fluxes for day 2002092 (P=196) and day 2007284 (P=68)
      DATA SEEFAC/1.0844,1.0959,1.1176,1.1364,1.1837,1.1418,1.2151,
     >           1.2826/
      !--  Test to see if need to scale - see DATRD2 subroutine      
      !.. Using a linear fit
      IF(NINT(UVFAC(58)).LE.0) THEN
         F107AV=(F107+F107A)*0.5
         DO LSR=38,45
           I=LSR-37
           UVFAC(I)=1.0
             A=(SEEFAC(I)-1)/128.0
             B=1-A*68.0
             UVFAC(LSR)=A*F107AV+B
             IF(UVFAC(LSR).LT.0.8) UVFAC(LSR)=0.8
             !WRITE(6,'(I6,9F11.3)') LSR,UVFAC(LSR)
         ENDDO
      ENDIF
      RETURN
      END
C:::::::::::::::::::::::::: FACSR_OLD ::::::::::::::::::::::::::::::::
C........ The Schumann-Runge factors are scaled according to F10.7
C........ from Torr et al. GRL 1980 p6063
      SUBROUTINE FACSR_OLD(UVFAC,F107,F107A)
      DIMENSION UVFAC(59),SRFLUX(8),SRA(8),SRB(8)
      !............. Schumann-Runge scaling
      DATA SRFLUX/2.4,1.4,.63,.44,.33,.17,.12,.053/
      !...... first two SRA and SRB values out of order in Marsha's paper
      DATA SRA/25.5,20.7,13.2,11.6,11.3,7.86,7.68,4.56/
      DATA SRB/222.,129.,53.4,36.0,25.0,11.3,6.35,2.05/

      !----  Test to see if need to scale - see DATRD2 subroutine      
      IF(NINT(UVFAC(58)).LE.0) THEN

         DO 505 I=38,50
           LSR=I-37
           UVFAC(I)=1.0
           IF(LSR.LE.8)
     >     UVFAC(I)=(SRA(LSR)*1.0E7*F107+SRB(LSR)*1.0E9)/SRFLUX(LSR)
     >         /1.0E11
 505     CONTINUE
      ENDIF
      DO I=38,45
        WRITE(6,'(I6,9F11.3)') I,UVFAC(I)
      ENDDO     
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE ZATOLT(CHI,G,D,E,TIME)
      !...... This routine converts solar zenith angle to local time
      H=((COS(CHI)-SIN(G)*SIN(D))/(COS(G)*COS(D)))
      IF(ABS(H).LE.1.0) THEN
         TIME=((ACOS(H)-E)*57.296/15.0)+12.0
      ELSE IF(H.LT.1.0) THEN
         TIME=24.0
      ELSE
         TIME=12.0
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE LTTOZA(TIME,GLAT,DECL,CHI)
      !...... This routine converts solar apparent time to solar zenith angle
      HH=(TIME-12.)*15*3.14159/180.
      CHI=ACOS(COS(GLAT)*COS(DECL)*COS(HH)+SIN(GLAT)*SIN(DECL))
      RETURN
      END
C:::::::::::::::::::::::::::: F2PEAK :::::::::::::::::::::::::::::::
C-- This subroutine determines the F2 peak electron density (NMF2)
C-- and height (HMF2) by fitting a parabola to the 3 points near
C-- the peak D1, D2, and D3 are the densities at H1, H2, and H3
 
      SUBROUTINE F2PEAK(D1,D2,D3,H1,H2,H3,NMF2,HMF2)
      DOUBLE PRECISION D1,D2,D3,H1,H2,H3,NMF2,HMF2
C--    Determine the coefficients A, B, C of the equation
C--    DEN = A * H**2 + B * H + C
C
C--  The first part takes care of low L values where the highest density
C--  may be at the equatorial point
      IF(NINT(H1).EQ.NINT(H3)) THEN
          NMF2=D2
          HMF2=H2
          RETURN
      ENDIF
C-- determine NMF2 and HMF2 for good data
      ANUMER = (D1-D2)*(H2-H3) - (D2-D3)*(H1-H2)
      ADENOM = (H1*H1-H2*H2)*(H2-H3) - (H2*H2-H3*H3)*(H1-H2)
      A= ANUMER/ADENOM
      B= (D1-D2 - A*(H1*H1-H2*H2))/(H1-H2)
      C= D1 - A*H1*H1 - B*H1
C--     Now the peak height and density
      HMF2 = -0.5*B/A
      NMF2 = A*HMF2*HMF2 + B*HMF2 + C
      RETURN
      END       
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE GTOB(IEXP,BLOND,BLATD,GLATD,GLOND)
C........ Conversion from geographic to centered dipole coordinates
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !....... convert to radians
      GLAT=GLATD/57.29578
      GLON=GLOND/57.29578
      !..... rotate in longitude
      XG=COS(GLAT)*COS(GLON-PLON)
      YG=COS(GLAT)*SIN(GLON-PLON)
      ZG=SIN(GLAT)
      !...... rotate in latitude
      XM=XG*SIN(PLAT)-ZG*COS(PLAT)
      YM=YG
      ZM=XG*COS(PLAT)+ZG*SIN(PLAT)
      !...... geomagnetic latitude and longitude in degrees
      BLATD=ASIN(ZM)*57.29578
      BLOND=(ATAN2(YM,XM))*57.29578
      IF(BLOND.LT.-0.01) BLOND=BLOND+360
      RETURN
      END
C::::::::::::::::::::::::: AMBBLK ::::::::::::::::::
      BLOCK DATA AMBBLK
      COMMON/AMBCOM/SECSAV,NY,NDY,LY,SX,ID,DECL
      DATA SECSAV,NY,NDY,LY,  SX,ID,DECL
     >     /0.0  , 0, 0 ,0 , 0.0, 0, 0.0  /
      END
C::::::::::::::::::::::::::: SOLDEC :::::::::::::::::::::::::
C...... This routine calculates the solar declination (DELTA radians)
C...... and ephemeris transit time (ETRAN in secs). ETRAN is the time
C...... the sun is overhead in Greenwich
C...... REFERENCE - page 484 "Explanatory Supplement to the Astronomical
C...... Almanac" Kenneth Seidelmann, University Science Books, 20 Edgehill 
C...... Road, Mill Valley, CA 94941, 1992
C...... IDAY is yyyyddd (eg 1988015 = feb 15), UT is the universal time
C...... in hours (0-24) 
       SUBROUTINE SOLDEC(IDAY,UT,DELTA,ETRAN)
      !.... Find day of year (0-366), and year (1996), then Julian day
       INTEGER JD,DAYNUM
       REAL T,YEAR,UT,L,G,LAMBDA,EPSIL,E,GHA,DELTA,SD,ETRAN,DTOR
       DOUBLE PRECISION DJD,DUT
       DATA DTOR/57.29578/   ! degrees to radians
      !..... Recover year and date and make sure UT is less than 24
       DAYNUM=MOD(IDAY,1000)
       YEAR=(IDAY-DAYNUM)/1000
       JD=INT(365.25*(YEAR-1900)+DAYNUM)+2415020      ! Julian day
       DJD=JD                                ! extra precision needed         
       DUT=UT
       T=(DJD+DUT/24.0-2451545.0)/36525.0         ! # of centuries
      !..     WRITE(6,*) DAYNUM,YEAR,JD,UT
       L=AMOD(280.460+36000.770*T,360.0)              ! aberration
       G=AMOD(357.528+35999.050*T,360.0)              ! mean anomaly
      !...... LAMBDA= ecliptic longitude. DELTA=obliquity of the ecliptic.
       LAMBDA=AMOD(L+1.915*SIN(G/DTOR)+0.020*SIN(2.0*G/DTOR),360.0)
       EPSIL=23.4393-0.01300*T
      !...... Equation of time. Time difference between noon and overhead sun
       E=-1.915*SIN(G/DTOR)-0.020*SIN(2.0*G/DTOR)
     >   +2.466*SIN(2*LAMBDA/DTOR)-0.053*SIN(4*LAMBDA/DTOR)
       GHA=15*UT-180+E                    ! Greenwich hour angle
       DELTA=ASIN(SIN(EPSIL/DTOR)*SIN(LAMBDA/DTOR))  ! solar declination
       SD=0.267/(1-0.017*COS(G/DTOR))    ! 
       ETRAN=(12-E/15)*3600 ! Ephemeris transit time in secs for FLIP
       RETURN
       END

C::::::::::::::::::::::::::::::::: SET_NOON ::::::::::::::::
C...... This subroutine resets the UT so that a run begins at noon
          SUBROUTINE SET_NOON(SEC     !.. Universal time in secs
     >                        ,PCO    !.. L shell
     >                        ,BLON)  !.. magnetic longitude
          IMPLICIT NONE
          DOUBLE PRECISION PCO
          REAL SEC,BLON,STIME,BLAT,GLATD,GLON,HRS
          INTEGER IEXP
          STIME=12.0
          IF(ABS(SEC).LT.24) STIME=ABS(SEC)
      !.. geographic longitude first using BTOG
          BLAT=57.3*ACOS(1.0/SQRT(PCO))
          CALL BTOG(IEXP,BLON,BLAT,GLATD,GLON)
          HRS=(STIME-(GLON)/15)
          IF(HRS.GT.24) HRS=HRS-24
          IF(HRS.LT.0) HRS=HRS+24
          IF(HRS.GE.23.0) HRS=0.1    ! ensures a longer settling period   
          SEC=HRS*3600
          !...WRITE(6,'(9F7.2)') SEC/3600.,BLON,BLAT,GLATD,GLON
          END
C:::::::::::::::::::::::::::::::::::: APEX_LVALUE ::::::::::::::::::::::::::::::::::::::::::::
C....... This routine sets up the apex grid based on altitude and geographic lat and long 
C....... It returns the L-value and magnetic longitude. It loads key field parameters into
C....... the FIELD_LINE_APEX_PARAMS module
C....... P. Richards July 2021
      SUBROUTINE APEX_LVALUE(YYYYDDD,glat,glon,alt,
     >              LVALUE,BLONG) !.. Output: L-shell and magnetic longitude
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      implicit none
      integer YYYYDDD 
      real date,glat,glon,alt,hr,qdlat,alon
      real bmag,si,malat    !.. Magnetic field, sin(DIP), Magnetic latitude
      real B(3),bhat(3)     !.. 
      REAL LVALUE,BLONG

      !.. Load key field parameters into the FIELD_LINE_APEX_PARAMS module
      glatX=glat 
      altX=alt    !.. glatX,glonX, and altX needed for other Apex routines 
      glonX=glon
      IF(glonX.GE.180.0) glonX=glonX-360.0   !.. Apex uses -180 to 180
            
      !.. Use initial station latitude to indicate -ve or +ve      
      !.. sign of glat indicates hemisphere for FLIP_APEX_GRID routine
      SIGNLAT=glat 
      
      !.. Determine the fractional date for apex_setup
      date=int(YYYYDDD/1000)+real(MOD(YYYYDDD,1000))/365.0

      IF(date.GT.2023) WRITE(6,'(/A,/A,/A,/A,/A,/A,/A,/A,/A)') 
     > ' *************************************************************',
     > ' *** Update the Apex.f90 file in 2025, 2030, 2035, etc.',
     > ' *** Do a web search for the latest copy at NCAR, replace the',
     > ' *** version in the FLIP source files and recompile.',
     > ' *** Then replace the .EXE file in folder FLIPEXE',
     > ' *************************************************************'
      !.. First setup the global Apex grid from -180 to 180 degrees.
      call apex_setup(date,altmax)

      !.. Determine the L-value and magnetic longitude from the modified Apex 
      !.. coordinates, quasi-dipole coordinates, base vectors,
      !.. and other parameters by interpolation from precalculated arrays.   
      CALL APEX_COORDS(
     >    glatX,glonX,altX,hrX,bmag,B,bhat,si,qdlat,alon,malat)

      !.. malat is the magnetic latitude      
      LVALUE=1.0/cos(malat/57.2958)/cos(malat/57.2958)    !.. PGR
      BLONG=alon
      RETURN
      END
C:::::::::::::::::::::::::::::::: APEX_COORDS :::::::::::::::::::::::::::::::::::::::
C....... This routine sets up the apex grid and returns the magnetic coordinates for
C....... a given altitude and geographic coordinates
C....... P. Richards July 2021
      !..    glat:    Geographic (geodetic) latitude (deg)
      !..    glon:    Geographic (geodetic) longitude (deg)
      !..    alt:     Altitude (km)
      !..    hr:      Reference altitude (km)
      !..    b(3):    Magnetic field components (east, north, up), in nT    
      !..    bhat(3): Components (E, N, V) of unit vector in B direction
      !..    bmag:    Magnitude of magnetic field (nT)
      !..    si:      Sine of magnetic inclination 
      !..    alon:    Apex long = modified apex long=quasi-dipole long(deg)
      !..    xlatm:   Modified Apex latitude (deg)
      !..    vmp:     Magnetic potential (T.m)
      !..    xlatqd:  Quasi-dipole latitude (deg)
      !..    W, D, be3, sim, F: See Richmond [1995] 
      SUBROUTINE APEX_COORDS(glat,glon,alt,hr,bmag,B,bhat,si,
     >      qdlat,alon,malat)
      use apex                    !.. Routines for the Apex magnetic field model
      implicit none
      integer ier
      real glat,glon,alt,hr,qdlat,alon
      real bmag,si,malat,vmp,W,D,Be3,sim,F
      real B(3),bhat(3),d1(3),d2(3),d3(3),e1(3),e2(3),e3(3),f1(3),f2(3)
      real f3(3),g1(3),g2(3),g3(3)

      ! Compute Modified Apex coordinates, quasi-dipole coordinates,
      ! base vectors and other parameters by interpolation from
      ! precalculated arrays. parameter alon was changed to qdlon
      call apex_mall(glat,glon,alt,hr,B,bhat,bmag,si,alon,malat,vmp,W,
     +  D,Be3,sim,d1,d2,d3,e1,e2,e3,qdlat,F,f1,f2,f3,g1,g2,g3,ier)
      RETURN
      END
C:::::::::::::::::::::::::::::::: FLIP_APEX_GRID ::::::::::::::::::::::::::::    
C...... Determine Apex parameters at each FLIP grid point
C...... P. Richards July 2021
      SUBROUTINE FLIP_APEX_GRID(JMAX,   !.. Maximum # of grid points
     >                             Z)   !.. FLIP altitude grid
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      use FIELD_LINE_APEX_GRID   !.. IDIM,gdlat,gdlon,bmagF,qdlat,qdlon,DIPX,DECX,COSDIPX
      implicit none
      integer I,J,JMAX
      real bmag,B(3),bhat(3),si,alon,malat,alat,qdlati,altF(IDIM)
      DOUBLE PRECISION Z(IDIM)

      !.. Compute Modified Apex coordinates, quasi-dipole coordinates,
      !.. base vectors and other parameters by interpolation from
      !.. precalculated arrays from apex_setup. 
      !.. This call establishes malat and alon for APEX_GEOLOC
      CALL APEX_COORDS(glatX,glonX,altX,hrX,bmag,B,bhat,si,qdlati,
     >   alon,malat)

      !.. Calculate the parameters for the FLIP field line grid. Note that 
      !.. the Southern Hemisphere is done in reverse order.
      DO J=1,JMAX
        I=J
        IF(SIGNLAT.LT.0) I=JMAX+1-J  !.. Reverse for southern station
        altF(I)=REAL(Z(I))
        !.. First determine geographic coordinates in each hemisphere
        IF(J.LE.(JMAX+1)/2) THEN
           CALL APEX_GEOLOC(malat,alon,altF(I),hrX,gdlat(I),gdlon(I))
        ELSE
           CALL APEX_GEOLOC(-malat,alon,altF(I),hrX,gdlat(I),gdlon(I))
        ENDIF
        qdlon(I)=alon      !.. save magnetic longitude

        !.. Now determine the magnetic field properties
        CALL APEX_COORDS(gdlat(I),gdlon(I),altF(I),hrX,bmagF(I),B,bhat,
     >     si,qdlat(I),alon,alat)
        DECX(I)=ATAN(bhat(1)/bhat(2))*57.2958
        DIPX(I)=ASIN(si)
        COSDIPX(I)=COS(DIPX(I))
        !.. FLIP COSDIP is -ve in Southern Hemisphere
        IF(DIPX(I).LT.0) COSDIPX(I)=-COSDIPX(I)
        !.. COSDIP should be +ve at the equator
        IF(J.EQ.(JMAX+1)/2) COSDIPX(I)=ABS(COSDIPX(I))
        DIPX(I)=DIPX(I)*57.2958
      ENDDO 
      RETURN
      END
C:::::::::::::::::::::::::::::::::: APEX_GEOLOC ::::::::::::::::::::::::::::
C....... Calculate the geographic latitude and longitude of a point
C....... on this field line
C...... P. Richards July 2021
      SUBROUTINE APEX_GEOLOC(malat,alon,alt,hr,gdlat,gdlon)
      use apex                    !.. Routines for the Apex magnetic field model
      implicit none
      real alt,hr,gdlat,gdlon,alon,malat

      !.. Given Modified Apex coordinates determine geographic coordinates
      call apex_m2g(malat,alon,alt,hr,gdlat,gdlon)

      RETURN
      END
C:::::::::::::::::::::::::::::::: BTOGAX :::::::::::::::::::::::    
C...... This routine replaces the old BTOG. It simply looks up the 
C...... geographic latitude and longitude from the field line grid
C...... P. Richards July 2021
      SUBROUTINE BTOGAX(J,GLAT,GLON,DEC,DIP)
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      use FIELD_LINE_APEX_GRID    !.. IDIM,gdlat,gdlon,bmagF,qdlat,qdlon,DIPX,DECX,COSDIPX
      IMPLICIT NONE
      INTEGER J                   ! Field line grid index
      REAL GLAT,GLON,DEC,DIP      !.. Latitude and longitude for output
      GLAT=gdlat(J)
      GLON=gdlon(J)
      DEC=DECX(J)
      DIP=DIPX(J)
      RETURN
      END

C:::::::::::::::::::::::::::::::: UPDATE_FLIP_GRID :::::::::::::::::::::::
      SUBROUTINE  UPDATE_FLIP_GRID(JMAX,Z,COSDIP,BM,
     >                       GL,GLATD,GLOND,GR)
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      use FIELD_LINE_APEX_GRID  !.. IDIM,gdlat,gdlon,bmagF,qdlat,qdlon,DIPX,DECX,COSDIPX
      IMPLICIT NONE
      INTEGER JMAX,I,IUNIT
      REAL GLATD(IDIM),GLOND(IDIM)
      DOUBLE PRECISION Z(IDIM),COSDIP(IDIM),BM(IDIM),GL(IDIM),GR(2,IDIM)
      IUNIT=17
      DO I=1,JMAX
        IF(IUNIT.EQ.-17.AND.(I.EQ.1.OR.mod(I,10).EQ.0)) 
     >  WRITE(IUNIT,'(4X,A,A,A)') 
     >  'I      alt       GLATD     gdlat     GLOND     gdlon',
     >  '          COSDIP         DEC           mag_fld',
     >  '              B_lat                DIP              GRAVITY'
        IF(gdlon(I).LT.0.0) gdlon(I)=gdlon(I)+360.0
        IF(IUNIT.EQ.-17) WRITE(IUNIT,'(I6,F10.3,15(2X,F8.2))') 
     >   I,Z(I),GLATD(I),gdlat(I),GLOND(I),gdlon(I),COSDIP(I),
     >   COSDIPX(I),DECX(I),BM(I),bmagF(I)/100000.,
     >   GL(I)*57.2958,qdlat(I),ACOS(ABS(COSDIP(I)))*57.2958,
     >   DIPX(I),GR(2,I),
     >   -SIN(DIPX(I)/57.2958)*3.98E+10/((RE+Z(I))**2)
       GLATD(I)=gdlat(I)
       GLOND(I)=gdlon(I)
       COSDIP(I)=COSDIPX(I)
      !.. GL(I)=qdlat(I)/57.2958
       BM(I)=bmagF(I)/100000.
       GR(2,I)=-SIN(DIPX(I)/57.2958)*3.98E+10/((RE+Z(I))**2)
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::::::: LVTOGX ::::::::::::::::::::::::::::
C....... When geographic lat and long are not specified, this routine uses 
C....... L-value and magnetic longitude to determine magnetic latitude and 
C....... routine apex_q2g to determine geographic coordinates from magnetic 
C....... coordinates.
C....... P. Richards July 2021
      SUBROUTINE LVTOGX(YYYYDDD,alt,LVALUE,qdlon)
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      implicit none
      integer ier,YYYYDDD
      real date,alt,gdlat,gdlon,qdlon,BLAT,LVALUE

      !.. Determine the fractional date for apex_setup
      date=int(YYYYDDD/1000)+real(MOD(YYYYDDD,1000))/365.0

      IF(date.GT.2023) WRITE(6,'(/A,/A,/A,/A,/A,/A,/A,/A,/A)') 
     > ' *************************************************************',
     > ' *** Update the Apex.f90 file in 2025, 2030, 2035, etc.',
     > ' *** Web search for the latest copy at NCAR, replace the',
     > ' *** version in the FLIP source files and recompile.',
     > ' *** Then replace the .EXE file in folder FLIPEXE',
     > ' *** See bottom of *LOG*.txt FLIP output file for more info',
     > ' *************************************************************'

      !.. First setup the global Apex grid from -180 to 180 degrees.
      call apex_setup(date,altmax)

      !.. Determine magnetic latitude thru alt km from the LVALUE 
      BLAT=57.2958*ACOS(SQRT((alt+6370.0)/6370.0/LVALUE))
      !.. Get geographic lat and long
      call apex_q2g(BLAT,qdlon,alt,gdlat,gdlon,ier)
      if(gdlon.LT.0.0) gdlon=gdlon+360.0   !.. Apex uses -180 to 180
c      WRITE(6,*) 'alt      BLAT  qdlon  gdlat  gdlon   L'
c      WRITE(6,'(9F7.2)')  altX,BLAT,qdlon,gdlat,gdlon,LVALUE

      !.. Load key field parameters into the FIELD_LINE_APEX_PARAMS module
      glatX=gdlat
      altX=alt    !.. glatX,glonX, and altX needed for other Apex routines 
      glonX=gdlon
      IF(glonX.GE.180.0) glonX=glonX-360.0   !.. Apex uses -180 to 180

      RETURN
      END
C:::::::::::::::::::::::::::::::::: SETUP_APEX ::::::::::::::::::::::::::::
C....... This routine sets up the Global Apex grid for the FLIP PRINT phase only
C....... P. Richards July 2021
      SUBROUTINE SETUP_APEX(YYYYDDD)
      use apex                    !.. Routines for the Apex magnetic field model
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT
      implicit none
      integer ier,YYYYDDD
      real date,alt,gdlat,gdlon,qdlon,BLAT,LVALUE

      !.. Determine the fractional date for apex_setup
      date=int(YYYYDDD/1000)+real(MOD(YYYYDDD,1000))/365.0

      !.. Setup the global Apex grid from -180 to 180 degrees.
      call apex_setup(date,altmax)
         
      RETURN
      END
C:::::::::::::::::::::::::::::::::: LTOGX ::::::::::::::::::::::::::::
C....... Use apex_q2g to determine geographic coordinates from
C....... given quasi-dipole coordinates.
C...... P. Richards July 2021
      SUBROUTINE LTOGX(alt,LVALUE,qdlon,gdlat,gdlon,ier)
      use apex                    !.. Routines for the Apex magnetic field model
      implicit none
      integer ier
      real alt,gdlat,gdlon,qdlon,LVALUE,BLAT

      BLAT=57.2958*ACOS(SQRT(1.0/LVALUE))
      WRITE(6,*) '  BLAT=',BLAT
      call apex_q2g(BLAT,qdlon,alt,gdlat,gdlon,ier)
      if(gdlon.LT.0.0) gdlon=gdlon+360.0

      RETURN
      END
C:::::::::::::::::::::::::::::::::: APEX_SZA ::::::::::::::::::::::::::::
c Example 5:  Test related routines subsol, cossza, magloctm, mlt2alon.
c  Ensure that IGRF coefficients are set up for the correct date
c  (since cofrm may have been called by apex_mka for a different date).
C...... P. Richards July 2021
      SUBROUTINE APEX_SZA(YYYYDDD,UTSECS,glat,glon,qdlon,sza,mlt)
      use apex                    !.. Routines for the Apex magnetic field model
      implicit none
      integer iyr,iday,ihr,imn,YYYYDDD,UTSECS
      real glat,glon,qdlon
      real sec,sbsllat,sbsllon,csza,sza,mlt 
    
      !.. Find subsolar latitude, longitude for a specific time:
      iyr=YYYYDDD/1000                 !.. Determine the year YYYY
      iday=MOD(YYYYDDD,1000)           !.. Determine the day DDD
      ihr=UTSECS/3600                  !.. Determine UT hour
      imn=(UTSECS-3600*ihr)/60         !.. Determine UT minutes
      sec=real(UTSECS-3600*ihr-60*imn) !.. Determine UT seconds
 
      !.. Find solar zenith angle at glat,glon for this time:
      call cossza(glat,glon,sbsllat,sbsllon,csza)
      sza = acos(csza)*rtd

      !.. Find magnetic local time
      call magloctm(qdlon,sbsllat,sbsllon,colat,elon,mlt)

      RETURN
      end



